### This script will perform Multiple Linear Regression on Microclimate and Topographic Data
### Also, it will produce plots of predicted tmax versus observed tmax

#Install package maSigPro

library('leaps')

# First define working directory and read in climate data

in.dir='/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/'
out.dir = '/home1/99/jc152199/Microclimate Statistical Downscale/Plots/solaralldata/'
setwd(in.dir)
a = read.csv('microMACROtopoVEG.csv',header=T)
b=read.csv('solar_regress.csv',header=T)
c=read.csv('solar_regress_all_days_newtopos.csv', header=T)
d=read.csv(paste(in.dir,'solar_regress_all_days_newtopos_2.csv',sep=""), header=T)
z=read.csv(paste(in.dir,'MicroMacroMinMaxASCII.csv',sep=''), header=T)

# Define a Linear Model framework for the Stepwise Regression

lm1 = lm(micro_max ~ AWAP_max + EVI + solar + coastdist, data = a)
lm2 = lm(micro_max ~ AWAP_max, data = a)

lm3 = lm(micro_max ~ AWAP_max + fpcmean + fpcvar + solar + coastdist + roaddist, data = c)
lm4 = lm(micro_max ~ AWAP_max, data = c)

#Produce partial plots for linear model

pdf(file=paste(out.dir,"check.pdf", sep=""))
crPlots(model=lm1, terms=~.)
dev.off()

# Plot predicted tmax versus observed tmax

png(file=paste(out.dir,"ObsVsPred.png", sep=""))

plot(a$micro_max, lm1$fitted.values, xlab = "Observed Max", ylab = "Predicted Max")

lines(x=c(0,50), y=c(0,50), lwd=1.5, col=3)
lines(x=c(predict(lm1,interval="confidence")[,1],predict(lm1,interval="confidence")[,3]), y=c(predict(lm1, interval="confidence")[,2],predict(lm1, interval="confidence")[,3]), lwd=1.5, col=2)

dev.off()

# Plot AWAP Max vs. Micro Max

png(file=paste(out.dir,"AWAPvMicro.png",sep=""))

plot(a$micro_max, a$AWAP_max, xlab="Micro Max", ylab="AWAP Max")

dev.off()

# Plot histograms of data from 'a'

pdf(file=paste(out.dir,"hist.pdf",sep=""))

hist(as.vector(a$AWAP_max))
hist(as.vector(a$EVI))
hist(as.vector(a$micro_max))
hist(as.vector(a$coastdist))
hist(as.vector(a$solar))

dev.off()

# Fucking around producing a Q-Q plot

pdf(file=paste(out.dir,"QQ.pdf",sep=""))
abline(a=0,b=1)
qqplot(a$micro_max, a$AWAP_max, main="QQ")
abline(a=0,b=1)
qqplot(a$micro_max, a$EVI, main="QQ")
abline(a=0,b=1)
qqplot(a$micro_max, a$coastdist, main="QQ")
abline(a=0,b=1)
qqplot(a$micro_max, a$solar, main="QQ")
dev.off()

# Produce plots of independent vs dependent

pdf(file=paste(out.dir,"XvsY.pdf",sep=""))

plot(a$micro_max, a$AWAP_max, xlab ="Micro Max", ylab ="AWAP Max")
plot(a$micro_max, a$EVI, xlab = "Micro Max", ylab = "EVI")
plot(a$micro_max, a$coastdist, xlab = "Micro Max", ylab = "coastdist")
plot(a$micro_max, a$solar, xlab = "Micro Max", ylab = "solar")

dev.off()

# Fucking around with BSS regression, this actually works and produces a model

t.subsets = regsubsets(micro_max ~ AWAP_max + fpcmean + fpcvar + solar + coastdist + roaddist, data = c,nbest=1)

summary(t.subsets)

best.model = data.frame(nvar.in.model=1:nrow(summary(t.subsets)$which),bic=summary(t.subsets)$bic)

data.names = summary(t.subsets)$obj$xnames[2:length(summary(t.subsets)$obj$xnames)]

bic.lm <- lm(as.formula(paste("micro_max~",paste(data.names[summary(t.subsets)$which[order(summary(t.subsets)$bic)[1],][-1]],collapse="+"))),data=c)

summary(bic.lm)

# Fucking around with BSS regression, this actually works and produces a model

t.subsets = regsubsets(micro_max ~ AWAP_max + fpcmean + fpcvar + solar + coastdist + roaddist, data = d,nbest=1)

summary(t.subsets)

best.model = data.frame(nvar.in.model=1:nrow(summary(t.subsets)$which),bic=summary(t.subsets)$bic)

data.names = summary(t.subsets)$obj$xnames[2:length(summary(t.subsets)$obj$xnames)]

bic.lm <- lm(as.formula(paste("micro_max~",paste(data.names[summary(t.subsets)$which[order(summary(t.subsets)$bic)[1],][-1]],collapse="+"))),data=d)

summary(bic.lm)

# Subsetting data set to only rainforest sites

e = d[-grep('31',d$site),]
f = e[-grep('32',e$site),]


# Fucking around with BSS regression, this actually works and produces a model

t.subsets = regsubsets(micro_min ~ AWAP_min + fpcmean + fpcvar + solar + coastdist + roaddist, data = z,nbest=1)

summary(t.subsets)

best.model = data.frame(nvar.in.model=1:nrow(summary(t.subsets)$which),bic=summary(t.subsets)$bic)

data.names = summary(t.subsets)$obj$xnames[2:length(summary(t.subsets)$obj$xnames)]

bic.lm <- lm(as.formula(paste("micro_min~",paste(data.names[summary(t.subsets)$which[order(summary(t.subsets)$bic)[1],][-1]],collapse="+"))),data=z)

summary(bic.lm)

# BSS Regression for Mins

t.subsets = regsubsets(micro_max ~ AWAP_max + fpcmean + fpcvar + solar + coastdist + roaddist, data = d,nbest=1)

summary(t.subsets)

best.model = data.frame(nvar.in.model=1:nrow(summary(t.subsets)$which),bic=summary(t.subsets)$bic)

data.names = summary(t.subsets)$obj$xnames[2:length(summary(t.subsets)$obj$xnames)]

bic.lm <- lm(as.formula(paste("micro_max~",paste(data.names[summary(t.subsets)$which[order(summary(t.subsets)$bic)[1],][-1]],collapse="+"))),data=d)

summary(bic.lm)





